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FORCE DENSITY FUNCTION RELATIONSHIPS IN 2-D 
GRANULAR MEDIA 

ROBERT C. YOUNGQUIST, PHILIP T. METZGER, AND KELLY N. KILTS * 

Abstract. An integral transform relationship is developed to convert between two important 
probability density functions (distributions) used in the study of contact forces in granular physics. 
Developing this transform has now made it possible to compare and relate various theoretical ap- 
proaches with one another and with the experimental data despite the fact that one may predict 
the Cartesian probability density and another the force magnitude probability density. Also, the 
transforms identify which functional forms are relevant to describe the probability density observed 
in nature, and so the modified Bessel function of the second kind has been identified as the relevant 
form for the Cartesian probability density corresponding to exponential forms in the force magnitude 
distribution. Furthermore, it is shown that this transform pair supplies a sufficient mathematical 
framework to describe the evolution of the force magnitude distribution under shearing. Apart from 
the choice of several coefficients, whose evolution of values must be explained in the physics, this 
framework successfully reproduces the features of the distribution that are taken to be an indicator 
of jamming and unjamming in a granular packing. 
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. 1. Introduction. A central topic within modern granular physics research is 
the study of intergranular force probability densities [1, 2, 3, 4, 5, 6]. The goal being 
to develop theory — from simplest assumptions — predicting the force density functions 
seen in simulations and experimental measurements. However, this goal is complicated 
by the differing forms for the force density functions presented in the literature, two 
of which are often treated as fundamental. 

The first of these force density functions is dependent upon the magnitude and 
angle of the contact forces between grains. It predicts force chains, the onset of 
granular jamming [7], strain hardening [1], and fracture of the individual grains. It 
has several odd and unexplained features, such as a finite, but non-zero, probability 
at zero force followed by an increasing probability to a peak near the average value of 
force. It has attracted scientific attention because its exponential tail at high forces 
and power law at weak forces are reminiscent of the Maxwell- Bolt zman distribution of 
velocities from statistical mechanics. Yet, its overall form is unlike any of the known 
statistical mechanics distributions and there have been numerous attempts to derive 
it [8, 9, 10, 11, 12, 13]. 

The second type of distribution that has been treated as fundamental in the 
physics literature is a force density function dependent upon the Cartesian components 
of the contact force vectors. This type of distribution is appealing because it relates to 
the following conservation law: If external forces are ignored and if the arrangement 
of grains is static then the total Cartesian force perpendicular to a plane cutting 
through -the -medium is conserved for any translation of the plane. Consequently, 
granular thermodynamic theories have been developed based on this conservation of 
the total Cartesian forces [14, 15]. 

The problem addressed by this paper is to clearly define these two types of force 
density functions and then to establish relationships between them. At first glance this . 
seems straightforward, and it would be if in the literature these densities were handled 
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as functions of two dimensions for two-dimensional media and three dimensions for 
three-dimensional media, but this is usually not the case. Probability force distribu- 
tions are typically expressed using one of the following variables; the force magnitude, 
the force angle, or a selected Cartesian component. Because of this contraction to a 
single variable pertinent information is not available and converting between the two 
types of densities is not possible. Yet, this conversion is important. Analytical theo- 
ries usually favor one distribution or the other and empirical investigations typically 
collect only one type of distribution and without a well established conversion it is 
not possible to compare the competing theories and their data against one another. 

The paper is organized as follows. Random variables are chosen such that the 
two two-dimensional probability force densities can be defined and the probability 
theory relations between them shown. In particular, an integral relation is developed 
whereby the two-dimensional polar force probability density can be converted to a 
Cartesian force density. From probability theory alone this relationship cannot be 
inverted but the integral relation can be recast as a set of Fourier transforms. Doing 
this allows an inverse relation to be found such that if the one- dimensional Cartesian 
force density function is known for all rotations of the axes, then the two-dimensional 
polar force probability function can be found. This inversion allows, for the isotropic 
case, the two force distributions to be treated as a transform pair. The properties 
of this transform relationship are discussed and significant solutions provided. An 
example is provided along with data generated from a Monte Carlo process. Then, 
an example of a pair of force density functions for the anisotropic case is given and 
compared with published results. 

2. Two-Dimensional Force Probability Density Functions. Suppose a 
large number of grains are placed randomly into a two-dimensional container and 
that the edges of the container push the grains against each other. These grains are 
not idealized and may be compressible, have arbitrary shape, be attracted or repulsed 
by each other, or have frictional contacts. The only restrictions on the grains are that 
they be static and that at each grain-to-grain contact a force vector be identifiable. 
Thus, the development presented here is applicable to a wide range of granular me- 
dia and could be extended to other discrete, static media such as intertwined fibers, 
foams, glass, and emulsions. Also, the contact forces between the grains and walls 
may be included in the density functions provided below as long as identifiable force 
vectors exist and that both force vectors, the grain on the wall and the wall on the 
grain, are counted in the statistics. 

At each grain-to-grain or grain-to-wall contact there are two force vectors, of 
equal magnitude and opposite direction, according to Newton’s third law, as shown 
in Figure 2.1. Assign the random variable, F, where 0 < F < oo, to the magnitude 
of these force vectors and the angle, 9 , where 0 < 6 < 27T, to the angle between 
them and the x-axis. A two-dimensional “polar” force probability density can then 
be defined in terms of the random variables, F and 0, and expressed as Pf,b{F, 9), 
describing the probability of finding a granular contact force with magnitude, F, and 
angle, 9. An immediate attribute of this function is that 


(2.1) P Fie (F,0)=P F>e (F i e+ tt), 

a result of there being equal and opposite forces at each granular contact (it is implied 
iu the definition of Pf,o(F , 9) that the angle 9 is modulo 27 t). 
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Fig. 2.1. A sketch showing intergrain forces and the associated random variables. 


A one-dimensional force magnitude probability density function, jPjp(F), can be 
defined by integrating over all values of 0 as shown in equation (2.2) below. This 
equation also shows that the symmetry relation given by equation (2.1) allows the 6 
integration to occur over any contiguous arc of length i r radians. 


a27t rn r<f>+ ^r/2 

(2.2) P f (F) = / P F , 9 (F,e)d6 = 2 / P Ft e(F,0)d9 = 2 / P F ,e(F,9) AO. 

Jo Jo J 4>—7r/2 

By integrating over all values of the force magnitude variable, a one-dimensional force 
angle probability density function, P$(0) is defined by 


(2.3) • 


Pe(6) = f 

Jo 


Pf,o{F>9) 


d F. 


From equations (2.1) and (2.3) the equal and opposite force result, i.e. Pe{9) — 
P e (Q- f 71-) can be shown. For frictionless grains (often studied in the physics literature) 
the direction of the .force vector is normal to the contacting surfaces of the grains. In 
that special case the density function of equation (2.3) is identical to the distribution 
of contact angles, referred to in the literature as the fabric of the granular material 
[16]. 

Not unexpectedly, there exists an alternative pair of random variables, F x = 
Fcos(0) and F y = Fsin(0), that can be used to describe a Cartesian force probability 
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density function, Pf x ,f v (Fx* F y ), which describes the probability of finding a granular 
contact force with Cartesian components, F x and F yy where 0 < F Xi F y < oo. The 
polar and Cartesian force probability densities are related through their corresponding 
Jacobians yielding 


(2.4) PfAM) = P Fm ,F v (Fcoa(0),Fsm(0))F 

and 

(2-5) Pf x ,f v (F x , F y ) = P F ,e{(F 2 x + Ft) 1 ' 2 , arctan (F y /F x ))/(F 2 + F 2 ) 1 / 2 . 

In equation (2.5), care should be taken to ensure that the arctangent returns an angle 
in the proper quadrant. Using equations (2.1) and (2.4), a Cartesian force density 
analogue for the existence of equal and opposite forces is found, 


(2-6) Pf x ,f v (F X! F y ) = Pf x ,f v (~F X! -F y ). 

Also, by integrating over either of the two Cartesian random variables, a one- 
dimensional Cartesian force density can be defined. Without loss of generality, in- 
tegrating over F y yields Pf x (F x ), which, as shown in equation (2.7) below, can also 
be found by integrating the polar two-dimensional density function expression from 
equation (2.5), 


(2.7) 


/ OO 

Pf x ,f v (F x ,F v ) d F y 

-oo 

°° PfA(F x + F 2 ) 1 / 2 , arctan {FJF X )) 


-L 


V2+f$) 1/2 


dF y . 


Changing variables from F y to # via tan 0 = F y /F x yields the following integral rela- 
tionship between the polar force density function and the one-dimensional Cartesian 
force density: 


) sec# d 9 


O ^tt/2 /*2tt \ 

1 + Pf,b(F x seed, 6) i 

0 J 3n/2 ) 

r3ir/2 

(2.8) Pf x (F x )= / Pf,o(F x sec#,#) sec# d# if F x < 0. 
Jn/2 


F x > 0 


Equations (2.6) and (2.7) imply that only one of the integrals in (2.8) needs to be 
calculated. 

A more general form for the Cartesian density function can be defined that will 
prove useful in the analysis that follows. Figure 2.1 shows a coordinate system rotated 
by angle (j), yielding new random variables, F x > and F y *. Recalling that rotation does 
not stretch space, the Jacobian of the transformation between the random variables 
F x ',F y ' and F xy F y is unity so their respective probability density functions are equal, 
i e. Pf x ,,f v , (F X ',Fy') = p f x ,f v (F x (F x , ,F y '),Fy(F x ',Fy'))- Using this result, a one- 
dimensional Cartesian force density along the x' axis, Pf x , (fv)> can be defined as 
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/ oo 

p F X ,,F V , (F X ',Fy>) dFy- 
-oo 


(2.9) 


/: 


PfA( F x> + Fy'Y 1 2 > axctan (F y ,/F x .) + 4 >) 

. (F$+F*y /2 


dF„ 


similar to the result in equation (2.7). Now, changing variables from F y * to 9 via 
arctan(F y //F x /) + = 9 , yields a result similar to that shown in equation (2.8) except 

that the additional (j) angle appears in the integrand and in the limits of integration. 
The variable transformation yields 


PfA f *’) = 


j r/2 p4>-\- 2tt \ 

= ( / + / \Pp t e{Fx f sec(0 - <£), 9) sec(9 — cj>) d9 if F x > >0 

\J <f> J <f>+37r/2 J 

r<f>- f 37r/2 

= / PF,e(Px' sec (9 — 0), 6) sec {9 — (j>) d 0 if F x > < 0. 


( 2 . 10 ) 


Even though the density function Pp x , (F x *) is well defined by the above equations it 
suffers from two shortcomings. Firstly, it is an explicit function of the angle <p , and 
this should be reflected in the notation chosen. Secondly, if the angle <f) is allowed to 
range from zero to 2ix radians there is an ambiguity in the choice of how to represent 
the domain of this function. Specifically, choosing a rotation angle, and a random 
variable, F x > , is identical to rotating by cj) + tt radians and choosing a value for the 
random variable of —F x f. We chose to resolve the second issue by the allowing the 
angle (j> to range from zero to 2ir radians and defining a new random variable, j F^, which 
is equal to F x > but is always positive, i.e. 0 < F+ < oo. The first issue above can then 
be resolved by adopting the notation Pp 4> (F ( p,4>) for the density function associated 
with this new random variable. The single variable subscript indicates that this is a 
one-dimensional density function but the two-dimensional domain shows that it is an 
explicit function of the two variables, F^ and <f>. So the function P F<f> (F ( f >i (f)) has the 
same “polar” domain as the density function 9). It can be explicitly expressed 

using the first integral expression above — where F x > > 0 — and merging the integrals 
by using the 2ir periodic nature of the variable 9 


r<p+lv /2 

(2.11) P F AF^,<j>) - 2 / PF,e(F<t,sec(0 - (f>),0)sec(9 - <j>) dO, 

J<f>-n/2 

where a factor of 2 has been added for normalization. Also, note that the function 
Pf# (F<f>, <^) 3 m order to be a density function, must be normalized at each angle, < p . 
In other words, the total number of forces does not change with the choice of cj), so 
the integral of Pp^F^, 0) over all F$ must always equal 1. 

Equation (2.11) is a general integral relation allowing integration of the “polar” 
two-dimensional probability force density function to yield any desired Cartesian pro- 
jection density function. This is a useful relationship within granular physics research 
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in that it allows the calculation of a Cartesian density function from a force magni- 
tude density function, but it is not a surprising result. Once the proper definitions 
are made the derivation is straightforward. The more difficult result is to perform the 
inverse operation, namely to find the polar form from the Cartesian form, but this is 
not possible within the realm of probability theory because the function, Pp# (F<p , 0), 
is not a two-dimensional probability density function. Even so, the inverse operation 
can be performed as demonstrated in the next section. 

Finally, by using the variable transformation, O' = 6 -f 7r, it can be shown that it 
can be shown that 


(2.12) Pf* (F^^ + n) = Pp* (F 4 ,, 4 >) 

3. Fourier Transform Representation of Force Density Integrals. In this 
section it will be shown that equation (2.11) above can be represented as a set of 
Fourier transforms. Accomplishing this allows equation (2.11) to be inverted by uti- 
lizing the Fourier transform inversion properties. This approach is similar to the 
mathematics used in tomography [17], but the development presented here is distinct 
for two reasons. First, the symmetry relations of equations (2.1) and (2.3) provide 
simplification that does not occur in tomography. Second, the emphasis in tomogra- 
phy is to generate three-dimensional functions from a set of two-dimensional images, 
while in the present development the goal is to obtain two-, or in a future work, 
three-dimensional representations of the force density functions from one-dimensional 
Cartesian projections. 

As a result of the symmetries obtained above in equations (2.1) and (2.12) we 
will require only restricted forms of the Fourier transforms. For example, we require 
only the Fourier cosine transform instead of the full exponential form. The Fourier 
cosine transform of a function f(x ) is expressed as 

F c [f(x); u] = J f(x) cos(ux) dx, 

which has the same form as the inverse cosine transform [18]. The Fourier Transform 
representations are written out in detail in order to resolve notational and normal- 
ization variations that appear in the literature. In addition to this one-dimensional 
transform, we will require a two-dimensional Fourier transform. In Cartesian form 
this is written as 


1 /»C>o rOO 

(3.1) Ti D[f(x,y)\(u,v)) = — / f(x,y)exp(i(ux + vy))dxdy, 

but since the density functions are in polar form, this transform will need to be 
expressed in polar form also. Using the change of variables x = F cos(6), y = 
Fsin(<9), u = Gcos(^), v — Gsin(<£), the following form for the two-dimensional, 
polar, Fourier transform is found, 


F2D[f{F,0)\ (G,<£)] 

/'OO 7T 


1 /‘OO e-ZTt 

= — / / f(F, 6) exp (iFG(cos(8) cos (<f) + sin(0) sin (<p)))F d F d8 
2tt Jo Jo 
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2t r 

f(F, 0) exp ( iFG cos (8 — <p))F d F d0, 

but we will only be transforming functions where f(F t 0) — f(F,6 + 7r). Using this 
symmetry simplifies this to a two-dimensional form of the Fourier cosine transform 



-I rOO r2-IT 

(3.2) Kd,c[W ey, (G, 4>)} =2 -J o J o f( F ’ °) C0 < FG co < d - W F dF de - 

Having established the above notation, now consider the following lemma which 
is the key result of this paper. 

Lemma 1 . The projected force density function is given by 

(3.3) Pf, 4>) = 2\/2^ c [T2d,c[PfA f > e )/ F > (G, 4>)]\ F <t>) ■ 

Proof 1. Note that 


2\Z2kF c \Ta D,c[ F F,e{P, 0)/F\ (G, <f>)]; P '$\ 

= 2 V 2 vyjl J™ [E jfjf P F fi(F, e) cos (FG cos(<9 - <j > )) dF d6> 

O /»2tt /'OO / rOO \ 

= £y J p F9 (F,e)^J cos(FG cos(0 - <j>)) cos(GF$) dGj dFd6 

n oo 

PfAF , 0) (5(Fcos(0 - 0) - F+)) dFdO. 

_ 

Changing the argument of the Dirac delta function , 8{x), is done using 


cos(GF < p) dG 


5{F cos{6 -$)- F+) = 5{F - F^ sec {6 - </>))/ 1 cos(0 - <j>)\. 

Thus integration with respect to F is well-defined . The limits of integration with 
respect to 0 can be compacted and shifted so that cos is always positive , removing 
the need for the absolute value and introducing a factor of two . Then 


2\p2snT c [T 2 d,c[ f f,o(F, 0)/F\ (G, <£)];F^] 

p<f>-\-ir/2 

= 2 P F ,e(F 4 > sec((9 - <j>), 6) sec(0 - <j>) d0 

J <f) — Tv/2 

The benefit of this form is that equation (3.3) can be immediately inverted because 
each of the Fourier transforms is its own inverse, yielding the second key result of this 
paper 

Lemma 2. The 2-D polar probability force density may be expressed as a function 
of the projected Cartesian force density by 

F 

2V2^ 


(3.4) 


Pf,b(F,0) = 


FiD,c[ F c[P Fi , (F<p, <t>)\ (G, </>)]; (F, 0)]. 
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Equivalently, 


(3.5) PfAW) 

y F<t> (^V> 4 ) cos (F<f>G) cos (FG cos (<f> — 9))G d G d<f> d F$. 

This integral can be integrated immediately with respect to the variable G , but this 
leads to an awkward functional form that is not amenable to solving using the standard 
look-up tables. Instead, we choose to leave it in this form allowing the order of 
integration to be carried out on. a case by case basis. This integral can be further 
integrated with respect to either F or 0 to yield the single variable force density 
functions, demonstrating that knowledge of the projected Cartesian force density 
function can yield an explicit form for the force magnitude density. 

Furthermore, for a frictionless packing, in which the force angles are the same as 
the contact angles as discussed above, one may integrate out F from equation (3.6) 
according to Equation (2.3) to obtain the fabric of the packing. This result is striking 
because it begins with a knowledge of only Pp <p (F^ ) , <j>), which is a set of contact force 
distributions that on the surface appear to have no contact angle information. This is 
the first time that force distributions, alone, have been directly related to the fabric 
of a packing and this may be important to developing a theory of granular rheology 
in which changes in the fabric and the forces are coupled. Also, if the empirical 
studies can determine a simple ^-dependence for Pp tft (F < p i <fr), then it may be possible 
to discern the fabric of a frictionless packing simply by sampling Pp^ at only a few 
orientations of <f>, or maybe only along the principle stress axes. 

4. Force Density Integrals for Isotropic Material. In this section the two 
integral equations derived above, equations (2.11) and (3.6), are simplified for the 
isotropic case yielding a pair of integral transform equations. Some of the properties 
of this transform pair are presented and a list of useful solutions are shown. We start 
with the following definition: 

Definition 4.1. An isotropic medium is one in which Pp i $(F ) 6 ) = Pp(F)/( 27t). 
Thus, an isotropic medium implies a force density with no angular dependence. As an 
immediate consequence of this definition, equation (2.11) can be simplified to show 
that an isotropic material also has no (j) dependence on its Cartesian projected force 
density, 


m oo 

- 


(2tt)' 


(4.1) Pp (Ft,<j>) = ~ f Pf(F<i, sec(0)) sec(0) d<9. 

n J-v/2 

So, for the rest of this section we will use the notation Pp^F^) for the Cartesian 
projected force density, since the function is no longer dependent upon the angle <p. 

Equation (4.1) can be put into a more useful form by changing variables from 
the angle 0 back to the force magnitude using F = F# sec (9). This yields the integral 
equation 


Pr* 


2 /°° Pf(F) 

' * > * Jf V OF 2 - F*y/ 2 


dF. 


(4.2) 
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Equation (3.6) can be simplified as well. Since neither of the density functions have 
angular dependence, the (j) integration on the right can be performed yielding the 
result 


(4.3) P f (F) = F rrPFtiFt) cos(F 0 G) J 0 (FG)G d F# d G 

Jo Jo 

where Jq is the zero-order Bessel function [[18], eqn. 3.715.18], where we choose to 
perform the integration with respect to F# before <2. 

Equations (4.2) and (4.3) are inverse transforms relating the force magnitude 
probability density, Pp(F), to the projected Cartesian force probability density, Pp^ (F^)| 
for isotropic two-dimensional granular materials. Such materials are of current theo- 
retical and experimental interest and these equations apply to any such material as 
long as an identifiable force can be assigned to the grain-to-grain, and, if included, 
grain- to-wall contacts. Consequently, a wide variety of force distributions are expected 
to result from current research, making it beneficial to discuss some of the properties 
of this transform pair as well as to present some of the more significant solution pairs. 

In the discussion below we will use the symbol to link transform pairs and we will 
not always normalize the pairs. 

Equations (4.2) and (4.3) are clearly linear in that if Pp^(F < f > ) <-► Pf(F ) and if 
Rp <f> (F^) +-* Rp(F) are two sets of solutions then 

(4.4) aPp, (F+) + bRp^ (JF*) «-> aPp(F) + bR F (F) 

is a solution pair, where a and b are arbitrary constants. Another method for gener- 
ating new solutions is to note that if Pp^ (F<p) «-► Pf{F) is a solution pair, then 


(4.5) 


dP P .(F*) dP F (F) 

F *-dF^~ F -dF~ 


is also a solution pair as a result of properties of the transform. By combining this 
result with the linearity result it can be shown that 


(4.6) 


a 2 P n (F,) 2 d 2 P F (F) 

F * dF% dF* 


is also a solution pair and more generally that 


(4.7) 


^PfAF*) t : F n d n Pp(F) 
F * dF$ dF- 


is a solution pair. This result is very useful for obtaining sets of similar solutions when 
trying to fit experimental or simulation data. 

A set of Bessel function integral equations [[18], eqn. 6.592.10, 6.592.12-15] pro- 
vides a method for obtaining useful solution pairs. Let Z v represent any of the Bessel 
functions, first kind J v , second kind Y Ui third kind H v , or modified second kind K u . 
Then after appropriate changes of variable the referenced integral equations can be 
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Table 4.1 

Solution pairs for Cartesian and force magnitude functions in isotropic packings generated from 
equation (4-8). 


Pf*(F*) 

P f (F) 

1. 

(2a/n)K 0 {aFj,) 

a exp (— aF) 

2. 

(2a 2 /it)F < i > Ki(ocF^) 

a 2 Fex p (—aF) 

3. 

(-F 4 ,F(2aM^(Ko(aF lf ,)) 

<t> 

a n-ri F n exp (—aF) 

4. 

(2a 3 /«)F 2 K 2 (aF<,) 

a 2 F( 1 + aF) exp (—aF) 

5. 

-Yo(aFt) 

cos(qF) 

6. 

MaF*) 

sin(aF) 

7. 

a exp (—aFf) 

a 2 FK 0 (aF) 

8. 

a 2 F < f > exp (-aF^) 

a 2 F(aFKi(aF) - K 0 (aF)) 

9. 

a ri->rlpn exp 

a 2 (-Fr^(FK 0 (aF)) 

10. 

(2/t r) cos(aF^) 

aF Jo(aF) 

11. 

(2/7r)sin(aF^) 

aFYo(aF) 


put into the form of either equation (4.2) or (4.3) yielding the following transform 
pair 


, A f~2~ Z v ^xli{otF^) Z v (aF) 

(4 ' 8) V 7ra F - 1/2 ~ F"~' ' 

Using this result and the differential generation result of equation (4.7), the solution 
pairs listed in table 4.1 can be found. Not unexpectedly, since the transforms derived 
above are between polar and Cartesian representations, the solution pairs are often 
a Cartesian function, i.e. a sine, cosine, or exponential; and a polar function, i.e. a 
Bessel function. Probability densities are positive functions with finite total integrals, 
so the exponential and modified Bessel function pairs are especially useful. Since 
the transform relationships are linear, it is worthwhile to show the cosine and sine 
solutions (both as Pf(F) and Pf+{F$)) so that, if desired, the Fourier components of 
a solution can be considered. 

Other solution pairs to equations (4.2) and (4.3) exist and can be found in the 
standard integral tables, but many of them can not be normalized. Table 4.2 shows 
five normalized solution pairs involving Modified Bessel Functions and Gaussians, 
which may be useful in granular physics applications. 

To close this section a normalized pair of solutions is shown that have been fit to 
the results of a Monte Carlo granular force model [12]. In this case, the grains are 
round, hard, frictionless disks in an isotropic packing. Over 20,000 grain forces were 
calculated, and the Cartesian force density function and the force magnitude density 
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Table 4.2 

More solution pairs for Cartesian and force magnitude distributions in isotropic packings. 



P f (F) 

1. 

2^~p(-aF*) 

2 olF exp (— aF 2 y 

2. 

2aF^exp(-aFp 

2aF(2aF 2 — l)exp(— aF 2 ) 

3. 

2 y/^exp (~aF 2 /2)K 0 (aF 2 /2) 

2^Jlxfn exp (— aF 2 ) 

5. 

(2a/ir)Kg (a.F^/2) 

(2a/ir)K 0 (aF) 

6. 

(2a 2 F 4 ,/ir 2 )K 0 (aF <t ,/2)K l (aF^/2) 

(2oc 2 F/tr)Ki {ocF) 


function were found empirically. A fit to the Cartesian force density function was then 
found using a three term modified Bessel function summation as shown below. The 
solution pairs given above were then used to determine the force magnitude density 
function. The expansion coefficients can be chosen so that both distributions provide 
an excellent fit to the empirical data, as seen in Figure 4.1, thus demonstrating the 
success of the transformation developed in this paper. The plotted density functions 
are 

(4.9) Pp^Ft) = C [UF*K 2 gV) - 2 F+Ki + 2 K 0 (|f,)] 

and 

(4.10) P f (F) = C [(11tt/2)F 2 + (11 - tt)F + tt] exp (~|F) 
where the normalization constant C = 7 r 2 / (128 — 6tt 4- 2tt 2 ). 



Fig. 4.1. (Left) Graph of Cartesian distribution from Monte Carlo simulation with analytical 
fit. The inset shows behavior F x — 1. (Right) Graph of force magnitude distribution from Monte 
Carlo simulation with analytical fit. The inset shows behavior F — 2. 

It should be helpful to further granular research that the modified Bessel Function 
of the second kind has been identified as the naturally-occurring form for the Cartesian 
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distribution, corresponding to the exponential forms of the polar distributions. The 
two “knees” in the curve that are visible in the Cartesian distribution of Figure 4. 1 
seem to be indicated in the Cartesian distributions of Bagi [14], as well, although this 
identification has not been previously made. 

5. An Example of an Anisotropic Solution Pair. In most real world cases 
granular media are subjected to greater total compressional force along one axis 
than the other, for example, in a gravity dominated situation. Consequently, the 
anisotropic case is of interest, although it is significantly more complicated. Fur- 
thermore, the physics of jamming and unjamming have emerged as possibly the key 
concepts in granula media. The anisotropic case is relevant to this because shear 
stress (an aspect of anisotropy in the stress state) is one of the three ways to un- 
jam granular media as represented by the jamming phase diagram [19]. It has been 
proposed that the evolution of Pp(F ) — from having a peak as shown in figure 2.1 to 
being monotonically decreasing — serves as an indicator of unjamming [7]. Numerical 
simulations have indeed shown that this evolution occurs when the material is un- 
jammed through stress anisotropy [1]. Therefore, an explanation for the evolution of 
Pp(F) is central to the aims of granular physics. In this section a reasonable form is 
selected for the density function Pp^ (F^, <t>) for a sample anisotropic case and then the 
integrals evaluated to yield Pp (F) . This solution is discussed and compared against 
published data [1]. The purpose is to demonstrate that this mathematical framework 
is a sufficient framework to include the evolution of Pp(F). 

The reason we start with Ff 0 (F^, <f>) instead of Pp t $(F,0) in this demonstration 
is because the Cartesian distribution is the one associated with the force conservation 
law, and it is through that conservation law that the anisotropy is injected into the 
problem. It is well known that the normal components (diagonal elements) in the 
stress tensor scale according to a xx = a*F 6 cos (2^) as the coordinate system is rotated 
through angle <j> . Hence, the quantity of conserved force normal to the layers of a 
granular material will also scale according to this form when the layer is oriented 
at angle <j>. The values of a and b are determined by the forces applied along the 
principle stress axes of the system. Based on the successful fits presented at the end 
of -the previous section we choose to let Pp^ (F^, <j>) be represented as a sum of the first 
three modified Bessel functions, but where this explicit angular dependence is added. 

2 / * \ n— 1 

(5.1 )Pf,{F^4>) = X>„ (jFr b ) (a + bcos(2<f>)) n+1 F£K n ((a + bcos(2(f>))F^,). 

The parameter b determines the amount of variation in force with angle, equaling zero 
for the isotropic case and approaching a for extreme anisotropy. Thus the force density 
is shifted towards higher forces along the T-axis with b nonzero. The (a-f 6cos(2<£)) n+1 
factor has been included to normalize the distribution at every particular value of <j> 
as required from the discussion above. The (a + b)/(a — b) factor is conjectural, to 
increase the weighting of the lower order terms as the anisotropy increases. Doing 
this yields results that correspond to dynamic simulations as seen in Figure 5.1, but 
whose basis is unclear. The point being that choices can be made for the Cartesian 
form of the force density functions, which can then be converted to force magnitude 
or force angle density functions (i.e. fabric) and compared to published data. 

The force magnitude density function is found by using this form for Pp <f> (F ( f, } <j)) 
in equation (3.6) and integrating with respect to 0 (see equation (2.2)) yielding 



FORCE DENSITY RELATIONSHIPS IN 2-D GRANULAR MEDIA 


13 


Pf(F) = 


F 

(2tt)2 



(a + 6 cos(2 4>)) n+1 Fg 


K n ((a + b cos( 24 >))F < f 1 ) cos (F$G) cos (FG cos (<f> — 9))G d G d d <f> d9. 


The integration with respect to 9 can be performed immediately yielding 2ttJq(FG), 
and the integration with respect to F<p can be performed via [18, eqn. 6.699.12], Then 
the integration with respect to G can be performed via [18, eqn. 6.565.4] yielding the 
partial result 


Pf(F) = 

E ^ - yjl (^) (a + bcos(2<f>)) n+3 ' 2 K n _ 1/2 (F(a + 6cos(2<£))) d<f>. 

Using the identities for the half order modified Bessel functions the integrations with 
respect to <j) can be made yielding the result, 

P F (F) = | exp (-aF)(a 0 (^) (a/ 0 (6F) - bh(bF)) + ai ((a 2 + b 2 )FI 0 (bF) 

(5.2) -(6 + 2abF)h(bF)) + a 2 ( ) (F(a 2 + 26 2 + a 3 F + 3 ab 2 F)I 0 (bF) 

\ a -f- b ) 

-(36 + 5abF + 3a 2 6F 2 + 6 3 J F 2 )/ 1 (6F))) 

where the exponential dependence on the force is expected from the isotropic case, 
but the Bessel function dependence on the parameter 6 is novel ( I (rr) is the modified 
Bessel function of the first kind). Using the values ao = 7t 2 /2, a\ = — 7r, q >2 ~ H> 
a — 7t/2 figure 5.1 shows a plot of equation (5.3) for various degrees of anisotropy. 

For 6 = 0, the isotropic case, the plot is identical to that shown in Figure 4.1 
and equation (5.3) reduces to equation (4.10), but as 6 increases the shape of the 
curve changes, slowly moving towards a pure exponential. This is in agreement with 
published simulation data where the force magnitude density function evolves in a 
similar fashion with increasing anisotropy [1]. This demonstrates that the mathemat- 
ical framework can produce this evolution naturally; nothing more exotic than the 
relative weighting of the Bessel terms need be invoked to produce it. 

6. Summary and Conclusions. It is possible within the straightforward tech- 
niques of probability theory to convert from Pp^F, 6) to Pp <f> (F ( j > , </>). Unfortunately, 
those techniques cannot provide a conversion in the opposite direction. However, we 
may recognize that the conversion in the forward direction is equivalent to the compo- 
sition of Fourier cosine transforms of the function. Since these transforms have their 
own well-defined inverses, the conversion from PF^F^.cj)) to Pf,q(F, 6) can likewise 
be expressed. 

This inverse conversion is interesting for several reasons. First, it allows theoret- 
ical models that only predict Cartesian force distributions (such as the q model [15]), 
to be directly compared against the force magnitude distributions, which have been 
more important to granular physics. Second, the inverse conversion indicates a previ- 
ously unrecognized relationship between the Cartesian force component distributions 
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Pf(F) 



Fig. 5.1. Evolution of force magnitude distribution with increasing anisotropy. The curve with 
the lowest probability density for F = 0 is the 6 = 0 isotropic case. The other curves , in order, 
correspond to b = 0.3, 0.6, 0.9, and 1.57. 


Pp 4> (F < f > , <j>) and the fabric of the granular packing, which may be important in future 
theoretical developments. Third, for the special case of isotropic granular packings 
in which the (j) and 6 dependencies may be eliminated, the transform pair reduces 
to a simple form that can be solved for a wide range of functions. This indicates 
which functional forms for the Cartesian components correspond to particular func- 
tional forms for the force magnitudes. Since it is well-known that the distribution of 
the latter has an exponential tail, the corresponding form for Pp x (Fx) ought to be 
modified Bessel functions of the second kind. Expansion in a series of such functions 
(of increasing order) display two characteristic “knees” when graphed, and indeed 
it turns out that such knees are observed in the empirical Cartesian distributions. 
Thus, the natural form for Pp x {Fx) appears to have been identified, and this should 
provide insight into the physical mechanisms that produce the distributions. Fourth, 
treating these modified Bessel functions with increasing anisotropic stress naturally 
produces an evolution of Pf(F) that depends upon the choice of coefficients in the 
series expansion. Prior research has associated this evolution with the occurrence of 
jamming and unjamming in granular packings, and so the inverse transform indicates 
that jamming may be described as in increase in weighting of the zeroth-order modi- 
fied Bessel function. This insight should be helpful to explain the physics of jamming 
and unjamming, which are important concepts in granular physics. Finally, we have 
noted that this inverse transform proves some of the integrals listed in Gradshteyn and 
Ryzhik to be valid over parameter domains larger than those indicated by Gradshteyn 
and Ryzhik. Because this inverse conversion identifies these relationships and natural 
functional forms for granular force distributions, its derivation should be a helpful 
contribution in future research into the physics of granular jamming and unjamming. 

The authors gratefully acknowledge William Miles of Stetson University for review 
and editing of this paper. 
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